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We study the fine-scale statistics of temperature and its derivatives in turbulent Rayleigh- 
Benard convection. Direct numerical simulations are carried out in a cylindrical cell with 
unit aspect ratio filled with a fluid with Prandtl number equal to 0.7 for Rayleigh num- 
bers between 10^ and 10^. The probability density function of the temperature or its 
fluctuations is found to be always non-Gaussian. The asymmetry and strength of devia- 
tions from the Gaussian distribution are quantified as a function of the cell height. The 
deviations of the temperature fluctuations from the local isotropy, as measured by the 
skewness of the vertical derivative of the temperature fluctuations, decrease in the bulk, 
but increase in the thermal boundary layer for growing Rayleigh number, respectively. 
Similar to the passive scalar mixing, the probability density function of the thermal dissi- 
pation rate deviates significantly from a log-normal distribution. The distribution is fitted 
well by a stretched exponential form. The tails become more extended with increasing 
Rayleigh number which displays an increasing degree of small-scale intermittency of the 
thermal dissipation field for both the bulk and the thermal boundary layer. We find that 
the thermal dissipation rate due to the temperature fluctuations is not only dominant 
in the bulk of the convection cell, but also yields a significant contribution to the total 
thermal dissipation in the thermal boundary layer. This is in contrast to the ansatz used 
in scaling theories and can explain the differences in the scaling of the total thermal 
dissipation rate with respect to the Rayleigh number. 



1. Introduction 

Convective turbulence which is driven by the action of buoyancy forces appears in a 
wide range of geophysical and astrophysical systems as well as in numerous technological 
applications such as chip cooling devices or indoor ventilation systems (Kadanoff 2001). 
The main focus of convection experiments has been on the precise quantification of the 
global turbulent heat transport through the cell which is measured by the dimension- 
less Nusselt number Nu as a function of the applied outer temperature difference, the 
properties of the working fluid and the geometry (see e.g. Niemela & Sreenivasan 2003; 
Funfschilling et al. 2005). The three dependencies are quantified by dimensionless outer 
parameters, the Rayleigh number Ra, the Prandtl number Pr and the aspect ratio T. 
The numbers are defined as 

agATH^ v D 

Ra= , Pr=-, r = — , (1.1) 

where a is the thermal expansion coefficient, g the gravitational acceleration, AT the 
outer temperature difference, H the height of the cell and D its diameter. The Prandtl 
number compares the kinematic viscosity of the fluid p to the thermal diffusivity k of 
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the temperature field. The experiments were able to explore a large range of Rayleigh 
numbers in the so-called hard turbulence regime from Ra ^ \{f up to Ra ^ 10^^ in 
the case of liquid helium (Niemela ct al. 2000). However, most experiments can provide 
pointwise measurements of time series of the turbulent fields only (see e.g. du Puits et al. 
2007 and references therein). Only recently, two-dimensional cuts through the fiow field 
have been analysed by combining particle image velocimetry and shadowgraph techniques 
(Xi, Lam & Xia 2004). 

In contrast, direct numerical simulations (DNS) are currently unable to reach in the 
very high Rayleigh number regime, in particular for aspect ratios F > 1. The advantage of 
DNS is however the fully resolved spatial and temporal information on the turbulent fields 
and local mechanisms of heat transfer (Kerr 1996; Verzicco & Camussi 2003; Hartlep, 
Tilgner & Busse 2005; Shishkina & Wagner 2006; Shishkina & Wagner 2007). 

Directly related to the Nusselt number is the mean of the thermal dissipation rate 
which is given by (e.g. Grossmann & Lohse 2000) 



The thermal dissipation rate field itself measures the magnitude of the temperature 
gradient and is defined as 



Here, x ~ {x, y, z) is the three-dimensional position vector and Xi ^ x,y or z. Both 
equations, (|1.2p and (|1.3p . imply that the statistics of the fluctuating thermal dissipation 
field is connected with the local fluctuations of the conductive heat transfer, as given by 
the local currents, j(x, t) = — /c(Vr). The measurement of the thermal dissipation field - 
and thus of spatial derivatives- is experimentally challenging, especially in high- Rayleigh 
number turbulence. Experiments can usually provide well-resolved time derivatives (Bel- 
monte & Libchaber 1996) or temporal increments of the temperature field (Zhou & Xia 
2002). In case of temporal increments in convection, closed forms of the corresponding 
probability density functions (PDF) of stretched exponential type have been constructed 
successfully by Ching (1991, 1993). For other flows, such as axisymmetry jets, turbu- 
lent channel flows or homogeneous isotropic turbulence, PDFs of the velocity increments 
have been constructed as a superposition of Gaussian distributions (Gastaing et al. 1990) 
or a product of a Gaussian random variable and a scale-dependent random multiplier 
(Ghevillard et al. 2006). In Rayleigh-Benard convection, time derivatives however cannot 
be translated into spatial derivatives by a Taylor frozen-flow hypothesis as in pipe or 
channel flows. Furthermore, for turbulence in a closed vessel, the concept of homogeneity 
is limited to the cell centre only. He, Tong & Xia (2007) were recently able to measure 
four temperature signals close to each other simultaneously to reconstruct temperature 
gradients. Their analysis disentangled contributions to the total thermal dissipation com- 
ing from the bulk and boundary layers. Rayleigh numbers Ra ^ 10^ were attained in the 
experiments, but the resolution of the gradients remained limited to scales larger or equal 
to the thermal boundary layer thickness. Direct numerical simulations by Kerr (1996) 
and more recently by Shishkina & Wagner (2007) focused on the geometric properties of 
thermal plumes, the structures that carry the heat away from the bottom plate. 

A detailed statistical analysis of the spatial derivatives of the temperature and the 
thermal dissipation field in different regions of the convection cell is thus still missing. 
In particular, knowledge about the strength of the fiuctuations around the mean tem- 
perature gradient and their dependence on the Rayleigh number allows for a validation 
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of predictions by scaling theories, e.g. those by Grossmann & Lohsc (2000) . Such analy- 
sis is also interesting from the perspective of passive scalar mixing in turbulence where 
progress in the understanding of the mechanisms that cause intermittent fluctuations has 
been made recently (Shraiman & Siggia 2000). A first and open point is to understand 
the differences between passive and active scalars such as the temperature in convective 
turbulence. For the passive scalar case, it is known that larger amplitudes of the dissipa- 
tion fields are mostly concentrated on fine scales (Kushnir, Schumacher & Brandt 2006) 
and that their statistical study puts rather large resolution constraints on DNS. This was 
discussed by Schumacher, Sreenivasan & Yeung (2005) and Schumacher & Sreenivasan 
(2005). 

In the present work, we provide a detailed height-dependent statistical analysis of 
the temperature, its fluctuations and spatial derivatives. Odd-order moments of spatial 
derivatives along an imposed outer mean scalar gradient have been used successfully 
to quantify deviations from local isotropy in shear flows or for the mixing of scalars in 
turbulence (Warhaft 2002; Pumir 1996; Schumacher & Sreenivasan 2003). We adopt these 
ideas for the present active scalar case and discuss the dependence of these anisotropy 
measures on Ra and in different regions of the convection cell. Finally, we study the 
statistics of the thermal dissipation rate in the bulk and cilosc; to the bottom and top plates 
of the cell. The results are then related to findings from the passive scalar dissipation 
field. Our study is intended to build a bridge between the mixing of passive and active 
scalar fields by comparing the statistical properties for both cases. The outline of this 
paper is as follows. In the next section, we present the equations of motion and details of 
the numerical model. The third section discusses the statistics of the total temperature, 
temperature fluctuations and their gradients. The fourth chapter is on the statistics and 
Rayleigh number dependence of the thermal dissipation field. Concluding remarks are 
given at the end. 



2. Equations of motion and numerical model 

The equations for an incompressible three-dimensional Navier-Stokes fluid in the 
Boussinesq approximation are solved in combination with the advection-diffusion equa- 
tion for the temperature field. The system is given by 

— + {u-\7)u = -Vp+iy\/^u + agTe^, (2.1) 

V-u = 0, (2.2) 

dT 

— + {u-V)T = kV^T, (2.3) 

where p(x, t) is the pressure, u(x, t) the velocity field and r(x, t) the total temperature 
field. The latter can be decomposed into a mean profile {T)a and fluctuations 0, 

T(x,i) - (^)A(z) + ^^(x,^). (2.4) 

Two different types of statistical averages are used. The average (•)y is taken over the 
whole cell volume and a sequence of statistically independent snapshots. The average 
(•)^ is calculated over the circular plane at a fixed height z and a sequence of snapshots. 
We thus combine always spatial and temporal averages. All quantities are expressed in 
characteristic units. Length scales are normalized with respect to the height of the cell, H, 
velocities with respect to the free-fall velocity, U = y/agATH, time scales with respect 
to H/U. Temperatures are given by < f = (T - Tc)/AT < 1 where Tc is the (cold) 
temperature at the top plate. When measured in these units, the thermal diffusivity 
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No. Ra Pr V Nu Re Nr N, Aa A St 



1 


1 X 10^ 


0.7 


1 


16.7 


682 


193 


97 


128 


0.0236 


0.0104 


0.0300 


2 


5 X 10^ 


0.7 


1 


25.9 


1444 


257 


129 


160 


0.0139 


0.0080 


0.0193 


3 


1 X 10** 


0.7 


1 


31.7 


2036 


271 


151 


200 


0.0111 


0.0069 


0.0158 


4 


5 X 10** 


0.7 


1 


52.4 


4530 


301 


201 


256 


0.0065 


0.0058 


0.0095 


5 


1 X 10'' 


0.7 


1 


65.0 


6255 


361 


181 


310 


0.0053 


0.0051 


0.0077 



Table 1. List of the dimensionless parameters Ra, Pr and F for different runs. The Reynolds 
number is Re — '^''"""^ and the Nusselt number is given by (|2.7|l and (12. 8p . Urms is an average 
over the whole cell volume and a sequence of statistically independent snapshots. All three 
velocity components, Ur, and Uz, are included. The Reynolds number defined in this way 
follows Re = Q.337?a° '"'^. In addition, we list the corresponding grid resolutions in azimuthal, 
radial and vertical directions {N,i,, Nr, N^), the grid spacing due to Grotzbach (1983), Ac, and 
the maximum of the geometric mean of the grid spacings for the present simulations as given 
by (|2.6p . A. 5t is the thermal boundary layer thickness, which is determined as H/{2Nu). 



becomes n — 1/ \/ RaPr and the kinematic viscosity = ^ Pr / Ra. Our studies are 
conducted for Pr — 0.7 in a cylindrical container of aspect ratio F = 1. The top and 
bottom plates have no-slip boundary conditions, u = 0, at a fixed temperature. The side 
walls are adiabatic no-slip boundaries, i.e., u = and dT /dr — 0. 

The equations are discretized on a staggered mesh and solved by a second-order finite 
difference scheme (Verzicco & Orlandi 1996; Verzicco & Camussi 2003). The pressure 
field p is determined by a two-dimensional Poisson solver after applying one-dimensional 
fast Fourier transformations in the azimuthal direction. The time advancement is done by 
a third-order Runge-Kutta scheme. The grid spacings arc non-equidistant in the radial 
and axial directions. In the vertical direction, they correspond to Tschebycheff collocation 
points. Table 1 lists the parameters of the present simulations: Ra, Pr, F, the resulting 
Reynolds and Nusselt numbers and the corresponding grid resolutions. The table also 
shows the grid spacing due to the resolution criterion of Grotzbach (1983) (see also Kerr 
(1996) for a discussion) which is given by 



for Pr < 1. The maximum of the geometric mean of the grid spacings in our runs is also 
provided by 



Figure [T] shows the vertical mean temperature profiles taken in three different ways 
in the turbulent flow for a Rayleigh number Ra — 10*. The standard mean tempera- 
ture profile {T)a{z) is compared with two vertical mean temperature profiles which are 
conducted at fixed radial positions and averaged azimuthally, {T(rQ, (j)))^(z). It can be 
seen that the slopes of the temperature vary significantly. While the profile {T)a{z) has 
almost zero slope at z = 0.5, a destabilizing slope, d{T)^{z)/dz < 0, is observed close 
to the sidewalk Close to the centre line, the situation is reversed, we observe a stabiliz- 
ing slope, d{T)^{z)/dz > 0. Our finding agrees with a recent experiment by Brown & 
Ahlers (2007) and confirms that the assumption of the temperature drop AT across the 
thermal boundary layer is a simplification. The presence of a non-vanishing local mean 
temperature gradient can have an impact on the statistics of the small-scale temperature 




(2.5) 



A = 



max 

X— (0,r,js) 



^/a^a^aT. 



(2.6) 
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Figure 1. Mean temperature profiles at Ra — 10*. The standard profile {T)a{z) is compared 
with two vertical temperature profiles which are calculated at fixed radial positions and averaged 
azimuthally. One profile is taken at the sidewall (ro = R), the other profile is close to the centre 
line (ro = 2Ar). All mean profiles vary from to 1 along the horizontal axis. The interval 
[0.4, 0.6] is shown only. 




Figure 2. Mean square fluctuations {0^)a as a function of z/H for Rayleigh numbers 
Ra = 10^, 10* and 10^ at an aspect ratio F = 1. The inset magnifies the thermal boundary 
layer. The solid horizontal lines correspond to 5t = H/{2Nu). 

fluctuations, similar to the passive scalar case. This point will be discussed further in 
section 3. 

The thermal boundary layer is resolved with at least 14 lateral grid planes for all 
runs. This is indicated in Figl2] where we plot the vertical mean square profiles of the 
temperature fluctuations and highlight the (non-equidistant) vertical grid spacing by the 
symbols in the inset. The horizontal lines correspond to the thermal boundary layer 
thickness estimated by 6t — H/(2Nu), which is always close to the position of the 
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(T-T^)/AT (T-T^)/AT 

Figure 3. Probability density functions (PDF) of the dimensionless temperature T. (a) All data 
are taken for Ra = 10^ and from the whole plane at height z. Heights are z = 0.55t ip), z = St 
(□), z = 25t (*), z = 4St (A) and z = 0.5H (+). (b) Data are taken at the centre plane z = 0.5 
for different Rayleigh numbers: Ra = lO'' (o), Ra = 10** (□) and Ra = 10^ (*). 




maximum of 9'^^^{z) — {9'^)a{z). All runs were performed for at least 80 large-scale eddy 
turnover times in order to gather sufficient statistics. We also verified that the vertical 
Nussclt number profile given by 

Nu{z) = AT ' (2.7) 

K ^ 

is a constant with a standard deviation of about 1%. The mean value of the Nusselt 
number profile, 

Nu^^f Nu{z)dz, (2.8) 
H Jo 

follows the scahng law Nu = No x Ra^^ , with No = 0.175 and /3 = 0.283. 
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Figure 6. Probability density functions (PDF) of vertical derivative of the temperature, dT/dz, 
taken in the same planes as in Fig. [5^. All data are again for Ra = 10^. Heights are z = 0.55t 
(o), z^St (□), 2 = 2St (*), z = 45t (A) and z = 0.5H (+). 
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Figure 7. Skewness of the vertical temperature derivative as a function of the Rayleigh 
number Ra. (a) Data analysis is performed in the center of the cell for points in 
V = {r,6,z\QAH < z < 0.6-ff}. Data are displayed in double-logarithmic axes. In addition, 
the decay laws of the skewness with respect to Ra are shown which would be in line with a 
return to isotropy (dashed lines), (b) Analysis is conducted in the thermal boundary layer for 
points in V — {r, 9 , z\ < z < St}- Data are displayed in logarithmic-linear diagram. 



3. Temperature and its derivative statistics 

3.1. Total temperature and its fluctuations 

Figure [3] shows the probabiUty density functions (PDF) of the dimensionless temperature 
T at different heights from the bottom plate and as a function of the Rayleigh number, 
respectively. The findings are in agreement with Siggia (1994) and Kerr (1996). As ex- 
pected, the PDFs in Fig. [3^ become increasingly symmetric toward the centre of the 
convection cell. This height-dependent asymmetry can be measured alternatively by the 
skewness of the temperature fluctuations, 6, which is given by 

= (3-1) 

The graphs are shown in Fig. 2^. A small skewness around 6t implies that rising and 
falling plumes are observed at the edge of the bottom and top thermal boundary lay- 
ers. Simultaneously, the colder and hotter pockets of temperature are generated in the 
vicinity of the thermal plumes which detach from the thermal boundary layer. This is 
a consequence of the incompressibility of the fluid as investigated more detailed in the 
Lagrangian frame of reference by Schumacher (2008). Up to a height of z/TJ « 0.1, the 
skewness increases monotonically before it declines to zero at the cell centre plane. The 
plots show that the maximum value of skewness increases with increasing Rayleigh num- 
ber. A large positive value of the skewness close to the bottom plate means that rising 
plumes with T{z) > {T) a{z) are dominant in this region of the convection cell. This is 
exactly the region in which a cusp is forming at T = 0.5 in the PDF (see Fig. [3^). 

The distribution of the temperature in the centre plane {z = 0.5) of the cell as a 
function of the Rayleigh number is shown in Fig.[3]3. The profiles are clearly non-Gaussian 
for all three Rayleigh numbers. The support of the PDF becomes narrower with increasing 
Ra. The PDF for the fluctuations 9 has the same shape in the centre plane. We found 
that the PDFs of 6 for the two larger Rayleigh numbers collapse when the temperature 
fluctuation is normalized by the corresponding root-mean-square value, 0rms, at z = 0.5. 
The PDF for Ra = 10^ has slightly smaller tails. The magnitude of the deviations from 
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Gaussianity can be measured by the flatness which is given by 




(3.2) 



Figure [IJd shows the vertical profiles of the flatness for the three Rayleigh numbers. 
The non-Gaussianity, i.e., a flatness larger than 3, is present for all three Ra. This is in 
agreement with other studies such as by Heslot et al. (1987), Castaing et al. (1989) or 
Ching (1991). It can also be observed that the flatness value levels off in the cell centre for 
the two larger Rayleigh numbers. This could indicate that the transition to the so-called 
hard turbulence regime of convection has been completed. Dimotakis (2005) discussed 
a mixing transition for passive scalars at Re ~ 10^. Above this threshold, a weaker Re 
dependence was predicted. A similar behaviour could be the reason for our observation. 
All flatness profiles have a minimum right above the thermal boundary layer. The slight 
asymmetry of the profile at Ra = 10^ is due to limitations in the statistical analysis. 
Even a longer time advancement and 140 statistically independent snapshots were not 
sufficient to obtain a symmetric profile since this analysis of the higher-order moments 
is done plane by plane. 

As already seen, the non-Gaussian PDF of the total temperature and the temperature 
fluctuations in the cell centre is a robust feature in thermal convection. It has been re- 
ported in experiments starting with Heslot et al. (1987) and Castaing et al. (1989). An 
almost exponential shape was emphasized by Yakhot (1989) on the basis of a hierarchy of 
moment equations for the temperature fluctuations. These ideas have been extended by 
Ching (1993) later. Yakhot suggested that this form prevails even for moderate Rayleigh 
numbers (which implies moderate Reynolds numbers), whenever regions with u^O > 
are dominant. A correlation between the vertical velocity and temperature fluctuation is 
nothing else but a flngerprint of a local blob of heat transferred through the cell. Inspec- 
tion of our data yields indeed that larger regions with Uzd > Q are present everywhere 
in the cell, even in the centre where coherent plumes are absent. This is shown in Fig. [5] 
where we plot the vertical dependence of the probability P{uz9 > 0) for three different 
Rayleigh numbers. It can be clearly seen that the probability decreases with increasing 
Rayleigh number, but remains significantly larger than 0.5 for all three cases. The plot 
is another manifestation of the net transfer of heat through the cell. 

Interestingly, the passive scalar turns out to be more sensitive to a particular driving 
and the Reynolds number of the adverting flow. When a mean scalar gradient is absent, 
the statistics is close to Gaussian or sub-Gaussian such as in Mydlarski & Warhaft (1998) 
or Watanabe & Gotoh (2004). In case of a non- vanishing mean scalar gradient, Pumir et 
al. (1991) suggested an exponentially distributed passive scalar as a generic feature, based 
on a simple one-dimensional random advection model. This was confirmed experimentally 
for sufficiently high Reynolds numbers (GoUub et al. 1991; Jayesh & Warhaft 1991; 
Jayesh & Warhaft 1992). However, sub-Gaussian distributions of 9 have also been found 
for this case by Overholt & Pope (1996), Ferchichi & Tavoularis (2002) and Schumacher 
& Sreenivasan (2005). Experiments by Gylfasson & Warhaft (2004) detected that the 
statistics depends on the particular driving of the turbulent flow. While an active grid 
caused sub-Gaussian fluctuations, a fine static grid caused super-Gaussian passive scalar 
distributions. Furthermore the deviations from Gaussian distributions decreased with 
increasing downstream distance in their experiments. The two findings imply that the 
ratio of the integral scales to the system scale is important. This is a point which has to 
be addressed more detailed. 
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3.2. Vertical derivative of temperature fluctuations 

An important question in turbulent convection is if the fluctuations of the temperature 
field in the bulk can be considered as being locally isotropic. In order to shed light on 
this point we adopt an approach that has been used successfully for passive scalars with 
a mean scalar gradient (Warhaft 2000). As we have seen in section 2, the local mean 
temperature gradient can deviate from zero while the overall profile sums up to a mean 
temperature gradient which is almost zero (see Fig. 1). Following Brown & Ahlers (2007), 
we define a ratio S that relates a difference of azimuthally averaged mean temperatures 
at fixed radial distance rg to the total temperature drop. It is given by 

^ (T(ro))^(;^ = 3H/4) - {T{r^))^{z = H/4) _ 

Close to the side wah, we get E{R) = -0.12 for Ra = lO"^, E{R) = -0.07 for Ra = 10^ 
and S(_R) = —0.04 for Ra = 10^. This large slope is mainly due to the rising and falling 
plumes. At the centre line, it follows S(2Ar) = 0.015 for Ra = 10^, E{2Ar) = 0.009 
for Ra = 10^ and E{2Ar) = 0.006 for Ra — 10^. Although decreasing with increasing 
Rayleigh number for both cases, their magnitude is not negligible. We observe that local 
mean temperature gradients are present. This supports the idea to conduct an analysis 
similar to passive scalar turbulence. 

The statistics of the vertical derivative of the temperature, dT/dz, is shown in Fig. 
\6\ The PDFs are analysed at the same heights as those for the temperature in Fig. [3^. 
A sharp decrease of the PDF at 0.5St is present in the vicinity of dT/dz — 0, simply 
because positive vertical derivatives are very unprobable close to the bottom plate. With 
increasing distance from the plate, the support of the PDF decreases. A negative skewness 
of all distributions can be observed which implies preferential derivatives along the mean 
negative temperature gradient. Recall also that a turbulent field is thought to be perfectly 
locally isotropic when all odd-order derivative moments are exactly zero. Therefore, the 
derivative skewness with respect to the temperature fluctuations 9, defined as 

SoAz)-ffi^^ (3-4) 
\\az) Iv 

is presented here. It measures deviations from the local isotropy at the smaller scales of 
convective turbulence. Figure [7] shows that such deviations are indeed present in different 
regions of the cylindrical domain. With increasing Rayleigh number and thus increasing 
Reynolds number, the derivative skewness Sd^e{z) decreases in magnitude in the centre 
region (see Fig.[7K). Following the original idea on the return to local isotropy in a simple 
shear flow (Lumley 1967) and its adaption to the passive scalar case (Warhaft 2000), a 
return would require a rather rapid decay law of the skewness with respect to the Taylor 
microscale Reynolds number Rx, i.e. \So^e \ ~ This was not found in passive scalar 
case for Schmidt numbers around unity (Warhaft 2000). Given the scaling dependence of 
the Reynolds number with respect to the Rayleigh number, which will be discussed more 
detailed in section 4.1, and the relation Re ^ a return to local isotropy in convection 
would require a decay law which is given by 

\SaM-Ra~'^^\ (3.5) 

for a fixed Pr. Here C is the scaling exponent for Re ^ Ra'^ relation, which was found 
to vary between ( = 0.43 and 0.49 in experiments (see section 4.1). Figure [7^ indicates 
that our data follow such a decay law for the larger Rayleigh numbers. This is in line 
with the observation that the local mean gradients in the bulk decrease in magnitude 



Fine-scale statistics in convective turbulence 



11 



with increasing Ra. The mean temperature gradient still causes ramps and cliffs of the 
temperature and thus a non-vanishing derivative skewness, but its impact seems to be 
weaker compared to a passive scalar which lacks a return to isotropy. A larger range 
of Rayleigh numbers and the analysis of the hyperskewness would allow to draw a firm 
conclusion. 

The opposite trend with respect to the Rayleigh number can be observed in the thermal 
boundary layer where the skewness Sg^g{z) is increasing in magnitude with Rayleigh 
number Ra (see Fig. [7}d) . Figure [5] shows that particularly the PDFs in the boundary 
layer are strongly skewed to negative derivatives. This observation is due to the fact that 
the heat transfer is mostly contributed by the second term of p.7[) there and that the 
fluctuations about the mean are consequently dominated by the conductive contributions 
jz — ~Hi{dO / dz). We will come back to this point in the next section when discussing 
the height dependence of the thermal dissipation. 



4. Thermal dissipation rate 

4.1. Vertical profiles and Rayleigh number scaling 

Figure [8] represents the isosurface plots of a snapshot of T (top) and of et (bottom) at 
Ra = 10^. We recognize the ridges in the temperature isosurface which correspond to the 
plumes that are detached from the bottom and top boundary layers. Associated with the 
rising and descending plumes are larger amplitudes of the thermal dissipation field which 
can be observed in the bottom panel, similar to findings by Kerr (1996) or Shishkina & 
Wagner (2007). The pattern of ridges (or stems) in the top panel is recaptured almost 
one-to-one in the bottom panel. The isosurface plot in the bottom panel indicates also 
that the local maxima of the thermal dissipation field et are dominantly found close to 
both plates and not in the cell centre. 

In order to quantify this vertical dependence more precisely, we decompose the thermal 
dissipation rate field into contributions which result from the mean temperature profile 
and the temperature fluctuations. With (|1.3p . it follows 



eT(x, t) 



\ az I oz oz 



(4.1) 



= e^T){z) + 2K^^^^+eg{-^,t) (4.2) 

The mixed term vanishes after averaging over planes at fixed height z and thus 

{eT)A{z)^e(T){z) + {eg)A{z). (4.3) 

Figure [9] plots the vertical profiles of both contributions to the total thermal dissipation 
rate for Rayleigh numbers Ra = 10^ (left column) and 10^ (right column). All values 
of the thermal dissipation fields are given in units of the corresponding ensemble mean, 
Nu/y/RaPr (see also (|1.2p ). The contribution of the mean temperature profile, e^x): is 
dominant in the thermal boundary layer and decreases rapidly to almost zero toward the 
cell centre (see panels (a) and (b)). This mean profile contribution increases with increas- 
ing Rayleigh number since the scale across which the significant temperature variation 
appears, St, becomes smaller. Conversely, the thermal dissipation rate due to the tem- 
perature fluctuations is dominant in the bulk and exceeds the mean profile contribution 
by more than 3 orders of magnitude there (not shown). As displayed in panels (c) and 
(d) of Fig.[9]for both Rayleigh numbers, the total dissipation remains significantly below 
Nu/y/RaPr in the cell centre. 
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Figure 8. Snapshots of the turbulent convection at Ra = 10®. Top: Isosurfaces of the tem- 
perature T are taken at 0.7 (bottom plate) and 0.3 (top plate), respectively. Bottom: The 
corresponding thermal dissipation field er is shown at the level 10{eT)v- 
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Figure 9. Vertical profiles of the two contributions to the total thermal dissipation rate which 
follow from the mean temperature and the temperature fluctuations (see (|4.2|l '). Data are for 
Rayleigh number Ra = 10^ (a,c) and Ra = 10^ (b,d). The three different terms have been 
normalized by the ensemble average of the total dissipation rate {eT)v = Nu/V RaPr. In panels 
(c) and (d), respectively, the vertical profiles of the total thermal dissipation rate are replotted 
in order to demonstrate that they decrease significantly below the global ensemble averages. 

From (|1.2p and the scaling relation Nu ^ Ra^ , the Rayleigh- number dependence of 
the ensemble average {€t)v yields 

{eT)v Ra^^^ = Ra^ , (4.4) 

which gives an exponent 7 — —0.217 in our simulations. This global scaling can be refined 
for different subvolumes of the cylindrical domain. Well in the centre, we define the bulk 
dissipation as 

2 



{^t)bu 



Ik 



dT 

dxi 



for X = (x, y, 0.4 < z < 0. 



ulk 



The thermal dissipation in the boundary layer is given by 

2 . 



\^t)bl 



dT 
dxi 



for — {x,y, z < 6t) • 



(4.5) 



(4.6) 



BL 



In line with the definitions of averages in section 2, we combine again an average over 
many realizations of the turbulent fields and a volume average which is taken now over the 
subvolumes BL and Bulk as defined in (|4.5|1 and (|4.6p . respectively. Following Grossmann 
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& Lohse (2000) (see their eq. (2.15)), one can estimate the thermal dissipation in the 
thermal boundary layer by assuming that it is due to the mean temperature drop across 
5t only and that {eT)BL ~ ^(t)- This leads to 

The factor St/H stands for the fraction of the total volume which is occupied by the 
boundary layer. The last step in (|4.7p holds if k is a constant. For Rayleigh numbers 
between 5 x 10'' and 10^, our DNS give a scaling with respect to Rayleigh number of 

(eT)sL « 0.016 X i?a°-i^ (4.8) 

This exponent is significantly smaller than the predictions by the scaling theories, (3 — 2/7 
or 1/3. We conclude that (|4.7p is not capturing the significant contribution from {eg)BL 
which can also be seen in Fig. [9l In addition, we find that (ee) bl displays almost the 
same scaling with Rayleigh number as {£t)bl in our direct numerical simulations. 

In the bulk of the cylindrical cell, the thermal dissipation rate is due to temperature 
fluctuations only and can be estimated by (see eq. (2.12) in Grossmann & Lohse (2000)) 

{eT)Buik « (ee)B«;fc 4^-^ - n^—^RePr , (4.9) 

Further explanation requires the knowledge of Re{Ra, Pr). Corresponding to experi- 
ments, an anomalous scaling dependence of the Reynolds number on the Rayleigh num- 
ber is found, i.e. deviations from a scaling exponent of C = 0.5 (Grossmann & Lohse 2002; 
Brown et al. 2007a). Furthermore, the definition of the Reynolds number, especially the 
choice of the characteristic velocity, leads to additional variations in the exponents. In 
summary, the range of exponents is found to be between C, = 0.43 and 0.49. It follows 
thus from that 

{^T)Bulk ~ Ra'^ , (4.10) 
for a fixed Prandtl number, as in the present study, with an exponent C > 0. However, 
our numerical simulations lead to 

(eT)i3«ifc « 0.32 X ii-a-o-^s , (4.11) 

which comes close to the experimental findings by He et al. (2007) and exhibits the re- 
ported trend of Verzicco & Camussi (2003). It does not vary significantly if the bulk 
volume fraction is increased in the analysis. For example, an exponent of —0.34 in (14. lip 
was obtained in the region 0.2 < z < 0.8. The Rayleigh number range, for which this least 
square fit was taken, is between 5 x 10^ and 10^. We conclude that the scaling exponent, 
which is by the way rather robust with respect to changes of the volume fraction, has the 
opposite sign in contrast to the scaling based on the estimate by Grossmann & Lohse 
(2000). The reason for this discrepancy could be that the dimensionally correct ansatz 
U{AT)'^ /H may not be appropriate for the thermal dissipation rate. Similar to fiuid tur- 
bulence, the Rayleigh numbers (and therefore the Reynolds numbers) are still too small 
to reach the regime where the thermal dissipation rate given in units of U{AT)'^/H is a 
constant. This is known as the dissipation anomaly (Donzis, Sreenivasan & Yeung, 2005). 
Both {eT)Buik/{UrmsT^„is/H) and {eT)Buik/{Ur,nsAT'^/H) stih decrease with increasing 
Ra for our data record {Urms defined as described in the caption of table 1). 

4.2. Stretched exponential behaviour in the tails 

The mismatch of the Rayleigh number scaling with the predictions by Grossmann & 
Lohse (2000) suggests that strong fluctuations of the thermal dissipation were not taken 
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Figure 10. Tails of the PDFs of the thermal dissipation rates of the temperature fluctuations, 
ee, in the bulk. Planes with z G [4:St,H — 45t] are included. The tails have been fitted by 
stretched exponentials as given by (|4.12p with the fit coefficients and exponents as indicated in 
the figure. The inset magnifies the PDFs for the smallest dissipation rate amplitudes in log-log 
scale. The data are normalized with respect to the rms value calculated over the whole volume 
and a sequence of statistically independent snapshots. 



into account in the scaling. We discuss therefore the statistics of eg in this section. Figure 
[TU] shows the PDF of the dissipation rate of temperature fluctuations eg in the bulk. The 
analysis is conducted in planes between z g [4(5r, H — ASt]- In correspondence with the 
passive scalar case (Overholt & Pope, 1996; Schumacher & Sreenivasan, 2005), we fit a 
stretched exponential to the fraction of the PDF which extends from the most probable 
amplitude (see the inset of Fig. [T0|) to the end of the tail. We define X = eg/{eg)v and 
X* = X — X„ip, with Xmp is the abscissa of the most probable (mp) value. Data are 
fitted to 

p(X*) = ^=exp(-mX*") . (4.12) 

Vx* 

Similar to the passive scalar case, we do observe fatter tails with increasing Rayleigh 
number (and therefore with increasing Reynolds number) . This is in line with an increas- 
ing degree of small-scale intermittency of the thermal dissipation rate field. Overholt & 
Pope (1996) and Schumacher & Sreenivasan (2005) have conducted a similar analysis for 
the passive scalar dissipation rate in a homogeneous isotropic turbulent flow. There, a 
minimal exponent a — 1/3 can be derived analytically in the limit of a very large Peclet 
number. The scalar is then advected in a flow which is white in time (Chertkov et al. 
1998, Gamba & Kolokolov 1999). The advection of a passive scalar in a Navier-Stokes 
flow always resulted in tails that were sparser, i.e. a > 1/3 (Yeung et al. 2005). Our 
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10 20 30 40 



V(^e)rms 

Figure 11. PDFs of ee at the cell centre (z=0.5H) normalized with respect to (es)rm3, which 
is also taken at z=0.5H. Present DNS data at Ra = 10^ are compared with the experimental 
resuhs (FIG. 4) of He et al. (2007) for Ra = 9.6 x 10** and Pr = 5.4. Statistics is gathered again 
over a sequence of independent snapshots. Data are more scattered in the far tail compared to 

Fig.iini 

present analysis indicates that this threshold is also approached for the hard-turbulence 
regime of convection. 

Figure fTD replots our data in the center plane for Ra — 10^ in units of the rms value 
of thermal dissipation {ee)rms m order to compare them with the measurements by He 
et al. (2007). A slightly sparser tail in the experiments is observed. On the one hand, the 
reason for this difference might be due to the improved resolution of the temperature 
gradient in our simulations, which was not possible in the experiment, as mentioned 
before. On the other hand, the Prandtl number {Pr = 5.4) was almost an order of 
magnitude larger in the experiment. An increasing Prandtl number reduces the thermal 
boundary layer thickness and can cause less plume mixing in the bulk. Less pronounced 
plumes will cause smaller amplitudes of the thermal dissipation. Such a trend with the 
Prandtl number would differ from passive scalar turbulence, where the tails of the PDF 
of ee become fatter with increasing Prandtl (or Schmidt) number at a given Reynolds 
number of the advecting flow (Schumacher & Sreenivasan, 2005). There a less diffusive 
scalar field generates sharper gradients. It remains open which of those two effects is more 
dominant. Further simulations are required to study the dependence of the temperature 
gradient statistics on the Prandtl number in thermal convection. 

When repeating the analysis for the part of the volume (z € [0,45t]) that has been 
excluded before, distributions as shown in Fig. [12] result. It is observed that the ampli- 
tudes of the dissipation increase by an order of magnitude for the largest Ra compared 
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Figure 12. Tails of the PDFs of eg close to the bottom plate {z e [0,45t]) for different Ra. 
The data are normalized with respect to the rms value calculated over the whole volume and a 
sequence of statistically independent snapshots. 



to Fig. [ini The fatter tails for the boundary layer data are in line with the analysis of the 
vertical plane averaged profiles (see Fig. The fluctuations of the thermal dissipation 
rate increase vi^ith increasing Ra. 

4.3. Deviations from log-normality 

The intermittent nature of the thermal dissipation field leads to deviations from log- 
normality, as can be seen in Fig. [13] where the same data set has been studied as in 
Fig. [lOl In order to highlight the differences in PDFs, we represent the data in linear- 
linear and linear-logarithmic axes respectively. It can be seen that neither the core nor 
both tails of the PDF fit perfectly. This is similar to the passive scalar mixing. Small- 
amplitude tails are fatter than the log-normal curve while large-amplitude tails remain 
sparser (Ferchichi & Tavoularis 2002, Schumacher & Sreenivasan 2005). The trend of 
the data with Rayleigh numbers between 10^ and 10^ suggests a very slow convergence 
towards the log-normality. The repetition of the analysis in the vicinity of the thermal 
boundary layers, i.e. for the same data set as in Fig. 1121 reveals even stronger deviations 
(see Fig. [14]). For the large-amplitude tail, a bump is observed. A similar feature has been 
reported by Kaczorowski & Wagner (2007). Our analysis suggests that this particular 
feature is due to the boundary layer dynamics. The bump in the distribution remains 
also when the PDF is analysed for the whole cell and not for a particular subvolume. 
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Figure 13. Replot of the same data as for Fig. [TU]in log-normal coordinates, fj, is the mean of 
In(ea) and o is the corresponding standard deviation. The solid line indicates a log- normal dis- 
tribution, (a) Linear- linear plot, (b) Linear- logarithmic plot of the same data. Rayleigh numbers 
are given in the legend. 
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Figure 14. Replot of the same data as for Fig. I12l in log-normal coordinates. Line and symbol 

styles are the same as in Fig. 1131 



4.4. Conditional mean thermal dissipation 

Deviations from the Gaussian distributions of the temperature fluctuation are also man- 
ifested in the conditional mean thermal dissipation. It is known that a conditional mean 
scalar dissipation is constant for a Gaussian scalar field, i.e. independent with respect to 
the scalar amplitude (see e.g. Overholt & Pope 1996). The mean is defined as 

p{ee,0{x,t) = ■0) 



{eo\ei^,t) ^ iP) = / degee ' 

Jo t) = V") 



(4.13) 



where p{es,9) is the joint PDF of the thermal dissipation rate and the temperature 
fluctuation. Figure [15] displays clearly the V-shape for the conditional mean for temper- 
ature fluctuation values around 0/9rms = 0. This is consistent with the deviations from 
Gaussianity, which were detected for the distributions of the temperature field. At larger 
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magnitudes of 9, the graphs reach a maximum before decreasing to zero since the support 
of the joint PDF and the temperature PDF is bounded. The V-shape is more pronounced 
with increasing Rayleigh number which is due to the higher thermal dissipation amph- 
tudes observable in the flow (see Fig.[TU]). Our results agree with the findings of Jayesh 
& Warhaft (1992). 



5. Conclusions 

We have presented a detailed numerical study of the statistics of the temperature fluc- 
tuations and their gradients in turbulent convection. The resolution constraints limited 
the study to moderate Rayleigh numbers between 10^ and 10^. We compared the fine- 
scale structure of the active temperature field to that of the passive scalar field, far away 
from the boundaries. In contrast to the passive scalar case, the temperature statistics 
is always non-Gaussian with a probability density function close to exponential in the 
bulk. Non-Gaussianity holds in all regions of the convection cell and for all mentioned 
Rayleigh numbers. Deviations in temperature fluctuations from the Gaussian distribution 
have also been confirmed by the V-shaped conditional mean dissipation. For the passive 
scalar fluctuations, the statistics depends on the mechanism that sustains the fluctuations 
and on the particular ratios of outer turbulence length scales. Super-Gaussian, Gausisan 
or weakly sub-Gaussian distributions have been observed there in the past. While the 
statistics of the temperature can differ from the passive scalar, the statistics of the spatial 
derivatives and dissipation rates behave qualitatively similarly. 
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The probability density function of the thermal dissipation rate in the bulk is fitted 
well with a stretched exponential. The tails extend to larger amplitudes with increasing 
Rayleigh number. This is a clear fingerprint of a stronger small-scale intcrniittcncy. Ther- 
mal dissipation field is always more intermittent in the boundary layer than the bulk of 
the convection cell. This is indicated by the extended tail of the PDF. All distributions 
showed clear deviations from the log-normality in the whole range of dissipation values. 

Motivated by similar studies in passive scalar turbulence and by our detection of locally 
varying mean temperature profiles, the deviations from locally isotropic temperature 
fluctuations are quantified by the third order vertical derivative moments as a function of 
the Rayleigh number. Interestingly, we found that the deviations from the local isotropy 
grow with increasing Ra in the boundary layer. In the bulk, the derivative skewness 
is found to decrease. Our data exhibit return to isotropy trend (Lumley 1967) when 
translating the original Reynolds number dependence of the skewness into a Rayleigh 
number dependence for the convection case. This is in contrast to the passive scalar case 
which does not show such trend (War haft 2000). The observed return might be due to the 
small amplitudes that the local mean temperature gradients possess. The vertical profiles 
of the plane averaged thermal dissipation show furthermore that an increasing fraction is 
concentrated in an ever thinner boundary layer with growing Ra. Finally, we demonstrate 
that the contribution of the temperature fluctuations to the total thermal dissipation in 
the thermal boundary layer is significant and cannot be neglected. This aspect is not 
included in the scaling of the thermal dissipation rate with respect to the Rayleigh 
number and can explain why the present trends of the total thermal dissipation with 
respect to the Rayleigh number differ in comparison to the scaling theory by Grossmann 
& Lohse (2000). 

Two aspects will be addressed in our future work on this subject. Firstly, an exten- 
sion to larger aspect ratios and higher Rayleigh numbers is desirable. It is known from 
experiments, such as by du Puits et at. (2007a), that the large-scale motion is sensitive 
to the aspect ratio. A rearrangement of the large-scale flow patterns will consequently 
affect the statistics of various quantities. Secondly, the direct link between the observed 
statistics and the local structures in the boundary layer and its vicinity is necessary. First 
steps have been done already in Shishkina & Wagner (2007) and Zhou et al. (2007). Our 
investigations on both aspects have started recently and will be discussed elsewhere. 
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